Load data

#> Info: preparing to read 1 data file(s)...
#> Info: reading file 1/1 'Garrett_Adducts.cf.rda' with '.rda' reader...
#> Info: loaded data for 71 data files from R Data Archive - checking loaded files for content consistency...
#> Info: encountered 6 problems in total.
#> # A tibble: 6 x 4
#>   file_id                  type  func           details                   
#>   <chr>                    <chr> <chr>          <chr>                     
#> 1 587__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
#> 2 587__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…
#> 3 586__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
#> 4 586__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…
#> 5 585__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
#> 6 585__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…

Vendor data table

#> Info: aggregating vendor data table without units from 71 data file(s), including file info 'c(ox = `Seed Oxidation`, everything())'

Map peaks

Chromatograms

Retrieve metadata

Map peaks

# add metadata
data_w_metadata <- vdt %>% 
  iso_add_metadata(metadata_samples, match_by = c(`Identifier 1`)) 
#> Info: metadata added to 1827 data rows, 40 left without metadata:
#> - 25 metadata entries were mapped to 1827 data entries via column 'Identifier 1'

# missing metadata
data_w_metadata %>% 
  iso_get_missing_metadata(select = c(Analysis, `Identifier 1`, map_id))
#> Info: fetching data entries that are missing metadata

# map peaks
data_w_peaks_all <- data_w_metadata %>% 
  # only continue with records that have metadata and are flagged for processing
  iso_remove_missing_metadata() %>% 
  filter(process == "yes") %>% 
  # map peaks
  iso_map_peaks(metadata_peak_maps, file_id = file_id, rt = Rt, rt_start = Start, rt_end = End)
#> Info: removing 40 of 1867 data entries because of missing metadata
#> 

# problematic peaks
prob <- data_w_peaks_all %>% 
  iso_get_problematic_peaks(select = c(file_id, Analysis, compound, Start, Rt, End))
#> Info: fetching 596 of 2044 data table entries with problematic peak identifications (unidentified, missing or ambiguous)

# focus on non problematic peaks only?
data_w_peaks <- data_w_peaks_all %>% iso_remove_problematic_peaks()
#> Info: removing 596 of 2044 data table entries because of problematic peak identifications (unidentified, missing or ambiguous)
data_w_peaks

Reference peaks

data_w_peaks %>% 
  iso_plot_ref_peaks(
    is_ref_condition = is_ref_peak == "yes", 
    ratio = c(`R 13C/12C`, `R 18O/16O`),
    group_id = `Analysis`,
    within_group = TRUE
  ) 

Data processing

Focus on analytes

Focus on the analytes and calculate a few summary parameters we want to use later.

Evaluate Standards

Known isotope values

Add known isotope values for standards.

#> Info: matching standards by 'type' and 'compound' - added 15 standard entries to 458 out of 762 rows

Single analysis calibration (for QA)

Determine calibrations fits for individual standard analyses.

#> Info: preparing data for calibration by grouping based on 'Analysis'
#> Info: generating 'd13C' calibration based on 1 model ('lm(`d 13C/12C` ~ true_d13C)') for 31 data group(s) with standards filter 'is_standard'. Storing residuals in new column 'd13C_resid'.
#> Info: there were no problematic calibrations
#> Warning in `[<-.data.frame`(`*tmp*`, is_list, value = structure(list(`3` =
#> "<>", : replacement element 1 has 1 row to replace 0 rows
#> Warning in `[<-.data.frame`(`*tmp*`, is_list, value = structure(list(`3` =
#> "<>", : replacement element 2 has 1 row to replace 0 rows

Calibration coefficients

Look at calibration parameters (coefficients and summary) as well ad data ranges in overview table.

Visualize the calibration parameters

# unnest some of the data for plotting
data_for_plots <- 
  stds_w_calibs %>%
  iso_unnest_data(
    c(Oxidation = `Seed Oxidation`, 
      datetime = file_datetime,
      `Mean Amplitude` = ampl_sample_mean.mV)
  ) 
data_for_plots


# parameters vs time
data_for_plots %>%
  iso_plot_calibration_parameters(
    calibration = "d13C", 
    x = datetime,
    color = `Injection Volume`,
    shape = Oxidation
  )


# parameter vs analysis
data_for_plots %>% 
  iso_plot_calibration_parameters(
    calibration = "d13C", 
    x = paste(Analysis),
    color = `Injection Volume`,
    shape = Oxidation
  )


# parameters vs amplitude
data_for_plots %>% 
  iso_plot_calibration_parameters(
    calibration = "d13C", 
    x = `Mean Amplitude`,
    color = `Injection Volume`
  ) 

# turn the last plot into an interactive one
ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

Inspect Residuals

Look at residuals to see goodness of fit for the single analysis calibrations.

stds_w_calibs %>% 
  # pull out all data
  iso_unnest_data(select = everything()) %>% 
  filter(is_standard) %>% 
  # calculate deviation from means
  group_by(Analysis) %>% 
  mutate(
    `Var: residual d13C [permil]` = d13C_resid,
    `Var: area diff from mean [%]` = (`Intensity 44`/mean(`Intensity 44`) - 1) * 100,
    `Var: amplitude diff from mean [%]` = (`Ampl 44`/mean(`Ampl 44`) - 1) * 100
  ) %>%
  ungroup() %>% 
  # visualize
  iso_plot_data(
    x = compound,
    y = starts_with("Var"),
    group = Analysis,
    color = `Injection Volume`
  ) +
  # plot modifications
  labs(x = "")

ggplotly()

Gobal calibration

Test different calibration models

Run global calibrations across all standards.

#> Info: preparing data for calibration by nesting the entire dataset
#> Info: generating 'd13C' calibration based on 4 models ('delta_only = lm(`d 13C/12C` ~ true_d13C)', 'delta_and_ampl = lm(`d 13C/12C` ~ true_d13C + `Ampl 44`)', 'delta_cross_ampl = lm(`d 13C/12C` ~ true_d13C * `Ampl 44`)', 'delta_and_time = lm(`d 13C/12C` ~ true_d13C + file_datetime)') for 1 data group(s) with standards filter 'type == "standard" & is_standard'. Storing residuals in new column 'd13C_resid'.

Visualize the calibration parameters for the different models.

# coefficient summary table
data_w_calibs %>% iso_unnest_calibration_parameters("d13C")

# plot to look at coefficients and summary
data_w_calibs %>% 
  iso_plot_calibration_parameters(
    calibration = "d13C",
    x = d13C_calib,
    color = d13C_calib,
    # include the significance level to gauge which model parameters matter
    shape = signif
  )

Garrett thinks delta and aplitude and delta cross amplitude look better

data_w_calibs %>% 
  iso_unnest_data(select = everything()) %>% 
  rename(`Injection Volume` = `Identifier 2`) %>% 
  filter(type == "standard" & is_standard) %>% 
  # plot each standard analysis for each of the models to see the differences
  iso_plot_data(
    x = compound,
    y = d13C_resid,
    group = paste(Analysis, d13C_calib),
    color = d13C_calib
  ) 

Apply calibration

Apply the delta_and_ampl calibration.

#> Info: applying 'd13C' calibration to infer 'true_d13C' for 1 data group(s); storing resulting value in new column 'true_d13C_pred' and estimated error in new column 'true_d13C_pred_se'. This may take a moment... finished.

Inspect results

The comparison of the predicted value in the standards along with the standard error in the prediction (the inherent error of the global calibration itself, which should be treated as the analytical error) can give a better sense of how precise the data is.

Visualization

Keep in mind that the calibration provides not just the predicted value of a peak but also the standard error based on the calibration and an indicator whether a predicted data point was in the calibration range or not.

data_calibrated %>% 
  iso_plot_data(
    x = Analysis, 
    y = true_d13C_pred, 
    y_error = true_d13C_pred_se,
    color = true_d13C_pred_in_range,
    shape = type,
    points = TRUE,
    lines = FALSE
  ) + 
  # additional features beyond the default plot
  facet_wrap(~compound, scales = "free_y") + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(color = "in calib. range")

Combine with other data for final plots

Remove bad chromtogram data

final_data <- data_calibrated %>%
  #filter out bad chromatograms - all are the lower amount injection, which were followed by a larger volume - keeping the larger volume injection data
  filter(Analysis != "597") %>% # 122 depth, OG005
  filter(Analysis != "593") %>% # 121.105, OG047
  filter(Analysis != "612") %>% # 121.89, OG040
  filter(Analysis != "613") %>% # 120.205, OG042-2
  filter(Analysis != "609") %>% # 119.55, OG158
  filter(Analysis != "595") %>% # 122.9, OG043
  filter(type == "sample") %>%
  filter(true_d13C_pred_in_range == TRUE) %>%
  select(
      # sample info
      depth = `Depth`,  `compound`, Identifier.1 = `Identifier 1`, `ampl_sample_mean.mV`,  type, Analysis, 
      # peak info
       Ampl.44 = `Ampl 44`,  true_d13C_pred , true_d13C_pred_se, true_d13C_pred_in_range,  d13C_resid
    ) %>%
  mutate(
    prep = "adduct"
  )

vis <- final_data %>% 
  filter(compound %in% c("C29", "C27", "C31", "C33"))%>%
  iso_plot_data(
    x = depth, 
    y = true_d13C_pred, 
    y_error = true_d13C_pred_se,
    size = Ampl.44
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse() +
  coord_flip()
ggplotly(vis)

Low res alkane data

lowres <- read.csv("lowres_d13c.csv") %>% select( -X, -`measurement_info`, -`file_datetime`, -file_id) %>% mutate(prep = "non-adduct") 


sapply(final_data, class)
#>                   depth                compound            Identifier.1 
#>               "numeric"             "character"             "character" 
#>     ampl_sample_mean.mV                    type                Analysis 
#>               "numeric"             "character"             "character" 
#>                 Ampl.44          true_d13C_pred       true_d13C_pred_se 
#>               "numeric"               "numeric"               "numeric" 
#> true_d13C_pred_in_range              d13C_resid                    prep 
#>               "logical"               "numeric"             "character"
sapply(lowres, class)
#>                compound            Identifier.1                   depth 
#>                "factor"                "factor"               "numeric" 
#>     ampl_sample_mean.mV                    type                 Ampl.44 
#>               "numeric"                "factor"               "numeric" 
#>               true_d13C          true_d13C_pred       true_d13C_pred_se 
#>               "numeric"               "numeric"               "numeric" 
#> true_d13C_pred_in_range              d13C_resid                    prep 
#>               "logical"               "numeric"             "character"

final_data <- bind_rows(lowres, final_data)
#> Warning in bind_rows_(x, .id): binding factor and character vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding character and factor vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding factor and character vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding character and factor vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding factor and character vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding character and factor vector,
#> coercing into character vector

Manipulate Alkane Data

Correct for biosynth fractionation

final_data <- final_data %>%
  filter(type == "sample") %>%
  select(-true_d13C, -d13C_resid) %>%
  filter(compound %in% c("phytane", "C17", "C19", "C21", "C27", "C29", "C31", "C33", "C35")) %>%
# net fractionation factor for inorganic -> lipid d13C (eps_net) calculated as = (inorganic -> biomass fractionation factor)  + (biomass -> lipid fractionation factor)

# minimum and max for marine algal biomass (first number in equation) from spread in data of algal biomass in Ohkouchi et al., 2015; OPTIONS FOR SHORTCHAIN: -4 per mil for phytane (Schouten et al., 1998; Sinninghe Damste et al., 2008), -8.4 per mil for cyano nC17 , -6.4 per mil for phytane (above takes into account frctionation in sulfur-bound and release per Sinninghe Damste?)(Sakata et al., 1997)
  
# lipid-leaf addition (second number in equation) by chain length from gymnosperm data in Diefendorf et al., 2011 [NOTE: average for gymnosperm is -2.5 per mil, that may be more appropriate]; minimum and maximum offset of leaf to atmosphere from C3 plant histogram in Ehleringer and Cerling, 2001 (in Taiz et al., 6th ed. Plant Physiology and Development)
  
  mutate(eps_net_min = case_when(.$compound %in% c("C17", "C19", "C21")  ~ (12+4) ,
                                 .$compound  == "phytane" ~ (12+4.4),
                                 .$compound  == "C27" ~ (22+2),  
                                 .$compound  == "C29" ~ (22+3),
                                 .$compound  == "C31" ~ (22+4),
                                 .$compound ==  "C33" ~ (22+3.5),
                                 .$compound ==  "C35" ~ (22+3)), 

         eps_net_max = case_when(.$compound %in% c("C17", "C19", "C21")  ~ (25+8) , 
                                 .$compound  == "phytane" ~ (25+6.4),
                                 .$compound  == "C27" ~ (32+2),  
                                 .$compound  == "C29" ~ (32+3),
                                 .$compound  == "C31" ~ (32+4), 
                                 .$compound ==  "C33" ~ (32+3.5),
                                 .$compound ==  "C35" ~ (32+3)),  
         
         d13C_pool_min = true_d13C_pred + eps_net_min,
         
         d13C_pool_max = true_d13C_pred + eps_net_max,
         
         d13C_pool_avg = true_d13C_pred + ((eps_net_max + eps_net_min) / 2) )

Calculate / estimate eps_p

final_data <- final_data %>%
    mutate(
         d13C_p = case_when(.$compound %in% c("C17", "C19", "C21")  ~ true_d13C_pred + 8.4,
                           .$compound == "phytane" ~ true_d13C_pred + 6.4 ) ,  # from Julio's unpublished text, original from Sakata et al., 1997 (6.4 in Sakata...)
         
         eps_p_lipid = 1000 * ((d13C_pool_avg + 1000)/ (d13C_p + 1000) -1), # from Freeman and Hayes 1992, but using avg estimated d13C of DIC from lipid rather than CaCO3 measurement. This allows moving pool value - important when pool itself is experiencing characterizable d13C changes thru time. 
         
         eps_f = case_when(.$compound %in% c("C17", "C19", "C21")  ~ 24,
                           .$compound == "phytane" ~ 25 ) ,    # from Julio's unpublished text, 22 for alkanes (cyano-plank intermediate) and 25 for phytane; Pagani 2002 says 25 for phytoplankton in general;  "Goericke et al. [ 1994] calculate the likely range of eps_f for marine phytoplankton to be 25-28" (from Brigidare et al., 1997). After Modelling (below), 24 for short chains gets us closest to phytane value (assuming 25 for phytane) - has implications for makeup of community, i.e., that cyanos are lower abundance than phytoplankton
         
         b = 170, # per Julio's unpublished text, can sequence from [118 to 262], or use constant 170. I think at least 170, potentially higher, given modelled need for high [PO4] advection to sustain OAE2
         
         CO2_aq_conc = b / (eps_f - eps_p_lipid) , # from Brigidare et al. 1997
         
         Ko = (2.418 * 10^-2), # solubility constant, from Weiss 1974 (tables II and III); dep. on T and Salinity, range from ~2 to ~7 over all regimes; @ ~35 C, range is from 2.1 to 2.3 across all salinities. Choosing for 34ppt (Hay et al., 1996) and 32C (mid range T from O'Brient et al., 2017)
         
         pCO2 = CO2_aq_conc / Ko  # from Julio's unpublished text , and Sarmiento and Gruber, Ocean Biogeochemical Dynamics, section 8.3 (p.327)
         )

Check calcultions immediately


final_data %>% filter(compound == "C21") %>% select( compound, eps_p_lipid, d13C_pool_avg, true_d13C_pred, CO2_aq_conc, pCO2) %>% kable()
compound eps_p_lipid d13C_pool_avg true_d13C_pred CO2_aq_conc pCO2
C21 16.45642 -5.558408 -30.05841 22.53572 931.9983
C21 16.44228 -4.716870 -29.21687 22.49355 930.2542
C21 16.44119 -4.652097 -29.15210 22.49031 930.1204
C21 16.41665 -3.188455 -27.68846 22.41754 927.1108
C21 16.42721 -3.818942 -28.31894 22.44880 928.4037
C21 16.42893 -3.921605 -28.42161 22.45390 928.6147
C21 16.39894 -2.129421 -26.62942 22.36531 924.9508
C21 16.47220 -6.495370 -30.99537 22.58295 933.9514
C21 16.45913 -5.719790 -30.21979 22.54383 932.3338
C21 16.45878 -5.698920 -30.19892 22.54278 932.2904
C21 16.43330 -4.182080 -28.68208 22.46687 929.1508
C21 16.45231 -5.313985 -29.81398 22.52345 931.4907
C21 16.39212 -1.721012 -26.22101 22.34527 924.1218
C21 16.42510 -3.692964 -28.19296 22.44254 928.1449
C21 16.39335 -1.794682 -26.29468 22.34888 924.2712
C21 16.40776 -2.656710 -27.15671 22.39127 926.0244
C21 16.43223 -4.118406 -28.61841 22.46370 929.0197
C21 16.48668 -7.354039 -31.85404 22.62648 935.7518
C21 16.44045 -4.608185 -29.10818 22.48812 930.0296
C21 16.45807 -5.656454 -30.15645 22.54065 932.2021
C21 16.44026 -4.596741 -29.09674 22.48755 930.0060
C21 16.44202 -4.701641 -29.20164 22.49279 930.2227
C21 16.44346 -4.787331 -29.28733 22.49707 930.3999
C21 16.43946 -4.548936 -29.04894 22.48516 929.9073
C21 16.38667 -1.394377 -25.89438 22.32927 923.4603
C21 16.43286 -4.155748 -28.65575 22.46555 929.0966
C21 16.45481 -5.462670 -29.96267 22.53091 931.7994
C21 16.45315 -5.363880 -29.86388 22.52595 931.5943
C21 16.44810 -5.063625 -29.56362 22.51090 930.9717
C21 16.44502 -4.880004 -29.38000 22.50170 930.5916
C21 16.44379 -4.806987 -29.30699 22.49805 930.4405
C21 16.43948 -4.550064 -29.05006 22.48521 929.9096
C21 16.45443 -5.439936 -29.93994 22.52977 931.7522
C21 16.49080 -7.598096 -32.09810 22.63890 936.2654
C21 16.45631 -5.551668 -30.05167 22.53538 931.9843

Model role of measured and fixed variables in pCO2 output

data_try <- expand.grid(
b = seq(170, 262, by = 10), 
Ko = seq(2.1, 2.9, by = 0.2) * 10^-2,
eps_p = seq(16.5, 17.5, by = 1)
) %>% as_data_frame() %>%
  mutate(
    eps_f = case_when(eps_p ==16.5 ~ 24,
                      eps_p == 17.5 ~ 25),
    compound = case_when (eps_p == 16.5 ~ "C17,19,21",
                          eps_p == 17.5 ~ "phytane"),
    CO2_aq = b / (eps_f - eps_p),
    pCO2 = CO2_aq / Ko
)

# base plot
base_plot <- 
  ggplot() +
  aes(x = b, y = pCO2, linetype = factor(Ko)) + #removed color for eps_p because both eps_p lie on the same trend (becuase of eps_f assignment)
  geom_line() 

  
# plotting data 1
base_plot %+% data_try

Bring in complimentary data (Jones et al., 2018)

Inorganic geochem data

cisotope <- read_excel(file.path("metadata", "Appendix_Table1_geochemistry.xlsx")) %>% 
  #rename columns
  select(depth = `Abs. depth (m)` , d13c_org = `Average δ13Corg (‰ VPDB)` )

final_data <- full_join(cisotope, final_data, by = "depth")

Fit to floating timescale


timescale <- read_excel(file.path("metadata", "Appendix_Table2_SH1_agemodel.xlsx")) %>% 
  rename(depth = `Depth`) %>%
  arrange(depth)

# extrapolation
extrapolated_timescale <-
  data_frame(
    depth = final_data$depth %>% unique()
  ) %>% 
  filter(!is.na(depth)) %>% 
  mutate(
    in_range = depth <= max(timescale$depth) & depth >= min(timescale$depth),
    Date.Ma = map_dbl(depth, function(d) {
      fit_data <- mutate(timescale, closest = abs(d - depth)) %>% arrange(closest)
      m <- lm(Date.Ma ~ depth, fit_data[1:10,])
      p <- predict(m, data.frame(depth = d), se.fit = TRUE)
      return(p$fit)
    }) ,
    in_range = depth <= max(timescale$depth) & depth >= min(timescale$depth),
    ka = map_dbl(depth, function(d) {
      fit_data <- mutate(timescale, closest = abs(d - depth)) %>% arrange(closest)
      m <- lm(ka ~ depth, fit_data[1:10,])
      p <- predict(m, data.frame(depth = d), se.fit = TRUE)
      return(p$fit)
    }))  
    
ggplot() + 
  aes(depth, Date.Ma) + 
  geom_line(data = timescale) +
  geom_point(data = extrapolated_timescale, mapping = aes(color = in_range))


alldata_w_time <- left_join(final_data, extrapolated_timescale, by = "depth")

Estimate error for deposition

alldata_w_time <-  alldata_w_time %>%
  mutate(age_range = case_when(.$compound %in% c("C17", "C19", "C21", "phytane")  ~ (-0.00001) , #in MA (100 yrs)
                                 .$compound %in% c("C27", "C29", "C31", "C33", "C35") ~ (-.01)), #in MA (10000 yrs)
         age_dev = Date.Ma + age_range 
   )

Final plots

High res


alkane <- alldata_w_time %>% 
  filter(compound %in% c( "C29", "C31", "C33"), true_d13C_pred_in_range == TRUE, prep == "adduct")%>%
  iso_plot_data(
    x = Date.Ma, 
    y = true_d13C_pred, 
    #y_error = true_d13C_pred_se ,
    shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse(limits = c(94.5, 94.34)) +
  #scale_x_continuous(limits = c(94.6, 94.3))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 
ggplotly(alkane)

bulkorg <-  subset(alldata_w_time, ka > -100 & ka < 120) %>%
  iso_plot_data(
    x = ka, 
    y = d13c_org
  ) + 
  # additional features beyond the default plot
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  #scale_y_reverse() +
  scale_x_continuous(limits = c(-100, 120)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()


grid.arrange(bulkorg, alkane, ncol = 2)

Low res (all)

alkaneall <- alldata_w_time %>% 
  filter(compound %in% c("pristane", "phytane", "C19", "C21", "C31", "C33"), true_d13C_pred_in_range == TRUE)%>%
  iso_plot_data(
    x = Date.Ma, 
    y = true_d13C_pred, 
    y_error = true_d13C_pred_se ,
    color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free", space = "free") + 
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse() +
  #scale_x_continuous()+
  scale_y_continuous() +
  coord_flip() 

ggplotly(alkaneall)

bulkorgall <-  subset(alldata_w_time) %>%
  iso_plot_data(
    x = ka, 
    y = d13c_org
  ) + 
  # additional features beyond the default plot
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  #scale_y_reverse() +
  #scale_x_continuous(limits = c(-100, 120)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()


grid.arrange(bulkorgall, alkaneall, ncol = 2)

Pool composition estimates

pools <- alldata_w_time %>% 
  filter(compound %in% c("C19", "C21", "C31", "C33"))%>%
  ggplot()+
  facet_grid(~compound, scales = "free", space = "free") + 
  geom_point(aes(x = Date.Ma , y = d13C_pool_avg)) +
  #geom_errorbar(aes(x = ka , ymin = d13C_pool_min, ymax = d13C_pool_max)) +
  #geom_errorbarh(aes(xmin = age_dev, xmax = ka , y = d13C_pool_avg)) +
  #geom_line(aes(x = age_dev, y = d13C_pool_avg, color = "red")) +
  #geom_line(aes(x = ka, y = d13C_pool_max)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "Inorganic pool d13C") +
  scale_x_reverse() +
  #scale_x_continuous()+
  scale_y_continuous() +
  #scale_x_continuous(limits = c(-100, 120)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()
ggplotly(pools)

Pool estimates with temporal offset estimates

pools_alttime <- alldata_w_time %>% 
  filter(Identifier.1 != "OG005_F1 alkanes_122m") %>% #removing duplicate for 94.45 Ma (non-adducted)
  filter(compound %in% c( "C19",  "C33"))%>%
  ggplot()+
  aes(x = Date.Ma , y = d13C_pool_avg, color = "red") +
  facet_grid(~compound, scales = "free", space = "free") + 
  geom_line() +
  #geom_errorbar(aes(x = ka , ymin = d13C_pool_min, ymax = d13C_pool_max)) +
  #geom_errorbarh(aes(xmin = age_dev, xmax = ka , y = d13C_pool_avg)) +
  #geom_line(aes(x = ka, y = d13C_pool_min)) +
  geom_line(aes(x = age_dev, y = d13C_pool_avg, color = "blue")) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "Inorganic pool d13C") +
  scale_x_reverse(limits = c(94.6, 94.3)) +
  #scale_x_continuous()+
  scale_y_continuous() +
  #scale_x_continuous(limits = c(94.6, 94.55)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()

ggplotly(pools_alttime)

CO2 calculations by date

alldata_w_time %>% filter(compound == "phytane") %>% select(Date.Ma, eps_p_lipid, d13C_pool_avg, true_d13C_pred, CO2_aq_conc, pCO2) %>% kable() 
Date.Ma eps_p_lipid d13C_pool_avg true_d13C_pred CO2_aq_conc pCO2
94.45489 17.93042 -6.505161 -30.40516 24.04670 994.4872
94.30791 17.92757 -6.349668 -30.24967 24.03699 994.0855
94.34293 17.91180 -5.490317 -29.39032 23.98351 991.8741
94.22115 17.88922 -4.257208 -28.15721 23.90736 988.7246
94.37841 17.87707 -3.592498 -27.49250 23.86659 987.0385
94.17981 17.90731 -5.245692 -29.14569 23.96835 991.2470
94.04433 17.87767 -3.625275 -27.52528 23.86860 987.1214
94.56738 17.91410 -5.615702 -29.51570 23.99130 992.1959
94.52951 17.91765 -5.809358 -29.70936 24.00333 992.6935
93.63407 17.92682 -6.309279 -30.20928 24.03447 993.9813
93.91610 17.90231 -4.972727 -28.87273 23.95147 990.5487
93.66267 17.92183 -6.037430 -29.93743 24.01752 993.2804
94.27372 17.89758 -4.714443 -28.61444 23.93552 989.8892
93.85049 17.89374 -4.504509 -28.40451 23.92258 989.3540
94.10260 17.88144 -3.831769 -27.73177 23.88124 987.6445
94.01944 17.87633 -3.551591 -27.45159 23.86409 986.9350
93.78417 17.90732 -5.246106 -29.14611 23.96838 991.2481
93.73495 17.92271 -6.085177 -29.98518 24.02049 993.4034
93.68296 17.91013 -5.399526 -29.29953 23.97788 991.6412
93.75582 17.91882 -5.872967 -29.77297 24.00728 992.8571
93.61487 17.91623 -5.731822 -29.63182 23.99851 992.4942
93.57814 17.90877 -5.324885 -29.22488 23.97326 991.4499
93.81802 17.91766 -5.810122 -29.71012 24.00338 992.6954
94.13923 17.89399 -4.517941 -28.41794 23.92341 989.3882

eps_p

#using my eps_p (lipd rather than caco3)
eps_p <- final_data %>% 
  filter(compound %in% c( "phytane" , "C19", "C21"), true_d13C_pred_in_range == TRUE)%>% 
  iso_plot_data(
    x = depth, 
    y = eps_p_lipid 
    #y_error = true_d13C_pred_se ,
    #shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "eps_p") +
  #scale_x_reverse() +
  #scale_x_continuous(limits = c(-100, 120))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 
ggplotly(eps_p)

pCO2

#using my eps_p (lipd rather than caco3)
pCO2 <- alldata_w_time %>% 
  filter(compound %in% c( "phytane" , "C19", "C21"), true_d13C_pred_in_range == TRUE)%>% 
  iso_plot_data(
    x = Date.Ma, 
    y = pCO2 
    #y_error = true_d13C_pred_se ,
    #shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "pCO2") +
  scale_x_reverse() +
  #scale_x_continuous(limits = c(-100, 120))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 
ggplotly(pCO2)

pCO22 <- alldata_w_time %>% 
  filter(compound %in% c( "C21"), true_d13C_pred_in_range == TRUE)%>% 
  ggplot()+
  aes(
    x = Date.Ma, 
    y = pCO2 
    #y_error = true_d13C_pred_se ,
    #shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_smooth() +
  geom_point() +
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "pCO2") +
  scale_x_reverse(limits = c(94.5, 94.34)) +
  #scale_x_continuous(limits = c(-100, 120))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 


alkane <- alldata_w_time %>% 
  filter(compound %in% c( "C31", "C33"), true_d13C_pred_in_range == TRUE, prep == "adduct")%>%
  ggplot()+
  aes(
    x = Date.Ma, 
    y = d13C_pool_avg, 
    #y_error = true_d13C_pred_se ,
    shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_smooth() +
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse(limits = c(94.5, 94.34)) +
  #scale_x_continuous(limits = c(94.6, 94.3))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 

ggplotly(pCO2)
ggplotly(alkane)
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
grid.arrange( alkane, pCO22, ncol = 2)
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).